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Steps for Solving the Radiative Transfer Equation for Arbitrary 
Flows in Stationary Spacetimes 
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ABSTRACT 

We derive the radiative transfer equation for arbitrary stationary relativistic flows in stationary 
spacetimes, i.e., for steady-state transfer problems. We show how the standard characteristics 
method of solution developed by Mihalas and used throughout the radiative transfer com- 
munity can be adapted to multi-dimensional applications with isotropic sources. Because the 
characteristics always coincide with geodesies and can always be specified by constants, di- 
rect integration of the characteristics derived from the transfer equation as commonly done 
in 1-D applications is not required. The characteristics are known for a specified metric from 
the geodesies. We give details in both flat and static spherically symmetric spacetimes. This 
work has direct application in 3-dimensional simulations of supernovae, gamma-ray bursts, 
and active galactic nuclei, as well as in modeling neutron star atmospheres. 

Key words: radiative transfer — relativity. 



1 INTRODUCTION 



The solution of the equation of radiative transfer in relativistic 
, flows is of considerable astrophysical interest. Steady-state and dy- 
' namical solutions of the transfer equation are particularly impor- 
tant for supernovae (SNe), gamma-ray bursts (GRB), and active 
galactic nuclei (AGN). The general f orm of the general relativis- 
tic transfer equation was derived by iLindguisl ( 1 1966) . who also 
derived the equation needed for neutrino transport in spherically 
synmietric flows such as the core-collapse o f massive stars . This 
work was further extended in Wilson ( 197lh. lBruennl ( Il985h . and 
iBaron et"al.i (1989). In the stellar community the fully-sp ecial rel- 
ativist ic transfer equation was derived and discussed by iMihafai 
( Il98(lh . General relativistic versi ons of the transfer equation hav e 
been derived and discussed by Morita & Ka nekol i 19841 1986h. 
[Schinder & Bludman ( 1989), Zane et al. ( 1996), as well bv lcastoil 
( l2004l) . lMihalasl (11980) proposed solving the transfer equation in 
the "comoving" frame using the method of characteristics. In this 
frame the momenta are measured by an observer moving with the 
flow, and the spatial coordinates are those of an inertial observer. 
In the steady-state spherically symmetric case the specific intensity 
is a function of only three variables, the inertial frame radius 
r, and two comoving momentum coordinates, the energy e = hv 
and fi the cosine of the angle between the direction of the photon's 
momentum and the radial direction. Comoving frames are partic- 
ularly useful because they simplify the form of the the collision 
terms in the transport equation. Mihalas made further use of the 
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spherical symmetry by treating the spatial and momentum angle 
variation of /,, separately from the dly/dv term. Somewhat confus- 
ingly he plotted "characteristic lines", which have one variable in 
real s pace and another in momentum space (see Fig. 1 of Mihala^ 
( ll980D ). These plots are curved lines and even though practitioners 
know that they are working in a mixed frame, it is common parlance 
to say "the characteristics are curved". Here we re-emphasize that 
photons move along geodesies (straight lines in flat spacetime) but 
demonstrate that for isotropic sources only one momentum vari- 
able (the energy) must be be comoving. We show that even in arbi- 
trary 3-dimensional flows one can choose parameters (coordinates) 
to label geodesies which do not change along phase-space char- 
acteristics (except for the afiine parameter, or "distance" along the 
characteristic itself). In addition we show that the change in comov- 
ing wavelength along the characteristic can be handled by stan- 
dard finite difference techniques. This procedure should simplify 
the development of fully 3-D radiation transfer codes both in flat 
space (applicable to variable stars, supernova and GRB spectra) 
and in curved spacetime (applicable to neutron star atmospheres 
and AGN). 

ISchinder & BludmanI ( Il989l) recognized that the momentum 
variables can be chosen as constants and the transfer equation sim- 
plified in the spherically symmetric case in the absence of a fluid 
flow. Although we developed our formulation independently, our 
work is an extension of theirs to the general 3-D case incorporating 
the effects of fluid flow. Their work used the method of variable 
Eddington factors, whereas our method is a characteristic based 
method and is specifically applicable to the case of arbitrary fluid 
flow. 
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In Sections |2] and [3] we introduce the Boltzmann and radia- 
tive transfer equations and the relevant phase-space quantities. In 
Sections ID and |5] we look at characteristics in flat and spherically 
symmetric spacetimes respectively. In Section|6]we discuss the log- 
ical steps necessary to solve the steady state transfer equation for 
stationary spacetimes and in Section|7]we give our concluding re- 
marks. 



2 THE BOLTZMANN EQUATION FOR PHOTONS 

In this section we do not restrict spacetime, but we neglect polar- 
ization effects. The Boltzmann Equation is an integro-differential 
equation for the invariant photon distribution function F{x, p) on 
the photon's 7-dimensional phase-space (x, p). This can be thought 
of as the photon"on-sheH" subspace of a full 8 dimensional particle 
phase-space. The number of photons AA' found by observer u{x) in 
a small 6-element AVjAP of phase-space at (x, p) is measured by 
the 6-form 6N, i.e., = 6N{hV„ AP) where 



6N = F(x,p)SV(„ 

and where (see lLindquislll 196^ for details) 
6V6 = -(u(x)- p)6V,,6P. 



(1) 



(2) 



In the above, u{x) is an arbitrary observer's unit 4- velocity at space- 
time point X, 



■ (ii{x) ■ p) = h/A 



(3) 



is the magnitude of the photon's 3-momentum as seen by observer 
u(x), 6 Vj is the observer dependent 3-dimensional volume element 
at X, and SP is the covariant volume element on the photons' 3- 
dimensional momentum space at (x, p). Here, h is Planck's constant 
and A is the wavelength measured by observer u(x). 

The coUisionless Boltzmann equation simply states that 
F[x(^), p(^)] remains invariant (constant) along the Lagrangian 
flow of photons in phase-space generated by their geodesic motion 
in spacetime. Constancy of F[x(^), p(^)] is a natural consequence 
of Liouville's theorem, i.e.,6V(, is invariant under this flow, and the 
constancy of AA' due to the absence of non-gravitational interac- 
tions of the photons. Any lack of constancy o f AA' in a finite vol- 
ume AVe is accounted for by a collision term jOxeniu not 
withstanding). To exhibit covariance, the Boltzmann equation with 
collisions, is often written as a differential equation 



dF 



dF\ 



_ dx" dF dp" dF 



(4) 



with F(x, p) explicitly given as a function of 8 variables (all 4 com- 
ponents of momentum are included but constrained by p • p = 0). 
The collision term on the R.H.S. is a measure of the rate of change 
of the number of photons AA' in a AVe transported along the would- 
be paths of non-interacting photons in phase-space. 

According to the geometrical optics approximation, photons 
travel on null spacetime geodesies independently of their wave- 
lengths. Affine parameters, ^, unique to each wavelength, can be 
chosen which generate the following orbits on phase-space: 



P" 



dx" 
dp" 

which reduces Q to 



(5) 
(6) 



' d.x" '^'^y'^P dp"-\d^)^, 



(7) 



The R.H.S. is typically separated into absorption and emission 
terms 



dF\ 



-fF + g, 



(8) 



where /(x, p) and g(x, p) are ide ntified respectively with the in- 
variant absorptivity and emissivity ("Morita & Kanekci|l986l) . These 
quantities implicitly depend on macroscopic properties of the in- 
teracting medium such as temperature, pressure, and density. The 
geodesic equations Jsj and ^ are equations for the characteris- 
tic curves of the integro-differential PDE Q. These characteris- 
tic curves (x(^), p(^)) are simply affinely parameterized spacetime 
geodesies, lifted to the 7-dimensional photon phase-space, which 
project back onto the null geodesies of spacetime x(^). It is impor- 
tant to understand that changing phase-space coordinates or chang- 
ing parameters for phase-space curves, doesn't alter these curves at 
all, only their description, e.g., when mixed coordinates are used as 
in Mihalas ( 1980) straight line geodesies naturally appear curved. 

Logically, solving the Boltzmann equation is a two step pro- 
cess: first solve the geodesic equations l|5j and l|6j for the required 
set of null geodesies and second solve the Boltzmann equation {JJl 
with appropriate boundary conditions. These two steps are com- 
bined in what is commonly called the characteristics method where 
Q is solved by changing the momentum variables to a comoving 
frame and the characteristic parameter to a nonphysical distance. 
These phase-space coordinate changes necessarily involve the fluid 
flow, and are unique only in highly symmetric cases. To tackle non- 
symmetric flows/spacetimes, however, the two steps are best kept 
separate. First solve the geodesic equations by using coordinates 
that are constant along the geodesies (or by finding the geodesies 
in any coordinate system and then transforming to constant coordi- 
nates) and second proceed to solve equation Q in the form of the 
transfer equation as given in the next section. 



3 THE RADIATION TRANSPORT EQUATION 

The radiation transport equation is an integro-differential equation 
for the specific intensity /(x, p) which is equivalent to the Boltz- 
mann equation l|4j or l|7j for F(x, p). Both are functions on the pho- 
ton's 7-dimensional phase-space (x, p); however, I(x,p) depends 
on a choice of an observer at each point of spacetime through^ 



I,{x,p) = —-{u{x)-pfF(x,p). 
h 



(9) 



We have chosen to follow our group's convention and use A rather 
than V as often appears in the literature; however, it is straightfor- 
ward to change between /( and ly using AI^ = vly. Once the ob- 
servers are chosen, Ia(x,p) like F(x, p), is a scalar. Defining a set 
of observers is equivalent to giving a unit time-like vector field on 
spacetime, m(x), which appears in Eq. Just as in equation (O, 



' /i(x, p)dAdAdCl is the rate observer u{x) detects energy crossing normal 
to his area dA in the direction p, within his solid angle dQ, and in his wave- 
band dA. Any locally-flat comoving reference frame at .v associated with the 
comoving observer u(x) can be used to evaluate dA and dil. These frames 
are arbitrary up to a rotation at each point x, however, actually defining one 
at every point x isn't necessary. 
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-u{x) ■ p is equal to the the photon's momentum as seen by ob- 
server u{x). If u{x) describes the material fluid with which the pho- 
tons interact, /( is called the comoving specific intensity and A the 
comoving wavelength. 

The transport equation for hix, p) is obtained from equations 
l|4j and l[8]l by substituting /^(x, p) for F(x, p) using Eq. ^ 



dh h 5dA h 



(10) 



where the observer dependent absorptivity xa and emissivity r] x are 
related to / and g by 

/ 



Xa 



(u(x) ■ p) ' 



(11) 



The other term (oc dA/d^) on the right in Eq. ilOl is present be- 
cause the definition of Eq. l|9), explicitly depends on comov- 
ing A. If as is customary we divide the extinction into two parts: 
"true absorption" and "scattering" cta, then^,i = + ctj. For 
a comoving observer we will also assume the emissivity is given 
by thermal emission (true absorption opacity ka times the Planck 
function B^) and that scattering is elastic and isotropic. This as- 
sumption is inherently required for our present formulation. It is 
beyond the scope of this work to consider anisotropic or inelastic 
scattering, but it is not entirely clear to us that the method cannot 
be extended to the more general case. For a comoving observer, xa 
depends only on the magnitude of the momentum and not its di- 
rection (given isotropic sources), and consequently is a function of 
only X and m • p in an arbitrary coordinate system. 

If the energy density in the radiation field is writte n as = 
AtUjJc (where J a is the classic 0* Eddington moment jMihalaj 
Il978l) . and Ea is defined in Eq. J47t). the emissivity becomes 



'7.1 



K,B, +cr,J, 



(12) 



When choosing coordinates (x, p) on phase-space for the purpose 
of evaluating I,i it is obvious that A itself should be one of the 
choices because it simplifies the evaluation of^^i. This is in fact 
the raison d'etre for using A as one of the momentum coordinates. 
When attempting a numerical solution, any dependence of xa on 
direction requires the use of a large number of angles in the kine- 
matically favored forward direction. Since as one moves around 
the interaction region, the forward direction changes, numerically 
resolving the variation of xa with direction can make the computa- 
tional requirements enormous. 

The reader should note that the right hand side of Eq. dlOl l 
differs slightly from the standard form of the non-relativistic static 
radiation transfer equation because the afflne parameter ^ is not a 
physical distance. As we discuss below it coincides with a distance 
(up to a constant) in some spaces, e.g., flat spacetime. 



4 FLAT SPACE SIMPLIFICATIONS 

When solving Eq. MQ\ for the comoving intensity /,i(x, p) or Eq. ([4]l 
for the Boltzmann distribution function F(x, p), the dimension of 
phase-space can effectively be reduced if there are common sym- 
metries in spacetime, the interacting medium, and the boundary 
conditions. This is because the transport equation has to be solved 
on only one of each equivalent set of characteristics on the full 
phase-space. For example if the spacetime, m(x), and the bound- 
ary conditions are stationary (i.e., have a timelike symmetry), a 



time dimension can be factored out and the required part of phase- 
space reduces to 6 dimensions (3 space and 3 momentum). If the 
spacetime, the flow, and the boundary conditions are stationary 
and spherically symmetric, phase-space reduces to 3 dimensions 
(1 space and 2 momentum). 



4.1 Arbitrary Stationary Flows 

For flat space with a stationary flow and static boundary condi- 
tions, 6 dimensions are required (3 space and 3 momentum). For 
space coordinates we choose the 3 Euclidean values r of the flat 
spacetime inertial system in which the flow is stationary 



u{x) = u(r) = r(r)(l,)S(r)). 



(13) 



For a boundary we choose a sphere of fixed radius R surround- 
ing the origin. We assume the emission and absorption coefficients 
vanish on and beyond r = R. On r = Rwe assume I,i(x,p) = for 
photons with incoming directions, i.e., with n • r < [see Eq. il4l 
below], and what we look for by integrating the transfer equation 
is the outgoing intensity on r = R, i.e., we seek Ix(x,p) for pho- 
tons with n • r > 0. To make Eq. MQ\ as easy to solve as possible, 
we choose two of the three momentum variables from the direction 
cosines of the photon's direction in the inertial system 



dx I dt dr\ h 



— (l,«\n\«'), 

A^ 



(14) 



e.g., we choose (n", n' ). The subscript "oo" refers to wavelength as 
seen by inertial observers. If we were to use A^ as a y' coordinate, 
all momentum coordinates would be constant, see Eq. l ll4b . and the 
characteristic equations l|5j and l|6j would be trivial. However, to 
accommodate the procedure used to solve Eq. l llOt we must use the 
comoving wavelength A as the third momentum coordinate: 



(15) 



(m(x)-p) 7(r)(l-n->S(r)) 

The comoving specific intensity Ia(x,p) = /,i(r, n) then de- 
pends on the spatial position r, two direction angles in n, and the 
comoving wavelength A; however, the comoving absorption and 
emission rj^ coefficients are independent of the direction angles. 
The obvious reason for choosing the first 5 of these 6 phase-space 
variables is that their dependence on the affine parameter ^ along 
any characteristic is easy to determine. The Euclidean position r(f) 
is linear in ^ and given by 



r(^) = b + i- n. 



where 



h 



(16) 



(17) 



The Euclidean distance .? is measured in the fixed inertial frame 
and b the impact vector defined by b • n = (see Fig. 1). If A„ is 
written as a combination of A and r using Eq. l ll5t . the combination 
in parenthesis in Eq. l ll6t still remains constant and the position co- 
ordinate part of a characteristic curve is still a straight line. The 
photon's direction n remains constant and although A(^) is a com- 
plicated function of f , it is given explicitly by substituting Eq. l ll6t 
into Eq. H5i . 

To logically connect with distant observations we can 
(somewhat arbitrarily) smoothly distort the comoving fluid u(x) to 
coincide with the rest observers at some point beyond r = R (see 
Fig. 1). Because and rj^ vanish in this domain, the fluid's effect 
on observations is sterile and A conveniently coincides with /loo. 




in out in out 



Figure 1. On the left a single fiat-space geodesic, Eq. U6t , described by the impact vector b and tangent n, enters and exits the boundary of the interaction 
region r = R. On the right, the related characteristic curves in phase-space con'esponding to a discrete set of wavelengths are also shown. If /t„ is used as the 
momentum coordinate the characteristics of Eq. (Ts) are all straight lines; however, if the comoving A is used the characteristics deviate from being straight 
but return to their A^o values beyond the boundary where the comoving fluid coincides with the rest observers. The characteristics of the differenced equation 
(43) are defined only when A„, is constant. 



In most applications, the transfer from comoving to stationary ob- 
servers is done abruptly at the boundary and is implemented by a 
Lorentz boost. The reason for using the comoving wavelength as 
the 6''' coordinate will become clearer when we alter the transport 
equation in such a way as to have characteristic curves which keep 
X constant. Rewriting Eq. l IlOl l by explicitly separating out the A 
dependence we obtain 



ldA\dIi I h 5dA\ h 



(18) 



If Euclidean coordinates on flat space are used as in Eq. (16} the 
first term in Eq. l |18l l is related to the r dependence of /,( by 



dh 



dr 

— • V/, 



(19) 



The dA/d^ part of the tangent to the characteristic is computed from 
Eq. ( list . If we would have chosen any 5 variables on phase-space, 
in addition to the comoving wavelength A, the transport equation 
would have been exactly of the form Eq. J18b . As long as the new 
coordinates for phase-space are given as functions of the origi- 
nal six, the characteristics in the new coordinates are found by 
direct substitution. The transfer equation in the form of Eq. dlSt 
says that /,i changes along a geodesic as usual because of emission- 
absorption and its explicit dependence on A due to its definition 
Eq. (|9]l; however, it also changes, when written as a function of A, 
because of its implicit dependence on a comoving wavelength that 
changes (A t 0) along the geodesic. Another way to say the same 
thing is that when geodesies on spacetime, are lifted to phase-space 
they are not constrained to A = constant hypersurfaces (see Fig. 1) 
and hence if A is used as one of the coordinates, Ia changes because 
A changes. 

4.2 Radial Stationary Flows 

In this section we try to clarify how our approach differs from the 
classic approach of Mihalas. To make contact with what is com- 
monly done in the stationary spherically symmetric case where the 



flow is radial = p{r)t and hence where phase-space reduces to 2 
independent dimensions beyond A, we introduce 1 inertial coordi- 
nate r and 1 comoving momentum coordinate /j. The wavelength 
and the radial coordinate can be evaluated from Eqs. dlSI l, ( |16| l, and 

GDI 

r(0 = iy- + s\ 

Vl -0-(r) 

AiO = /loo . (20) 



The coordinate yu is the cosine of the angle between the radial di- 
rection and the direction of the photon in a frame instantaneously 
moving with the radial flow at r and is found by a radial Lorentz 
boost, 



n-f-/?(r) sir-pir) 



(21) 



l-/?(r)n-f l-/3(r)i/r 

Equations l l20t and ( I21t are the integrated characteristic equation 
for the resulting transport equation for Ii{r,ji), 



drdh djjdlA dAdh_ I h. 5 dA\ h 



(22) 



which is equivalent to equation (3.1) of lMihalasI il98(t) for Mi- 
halas' parameter sm is related to the standard affine parameter <f by 
dsM = -{u- p)d^ = {h/A)d^ = (Ac^lX)ds, and equals the differential 
spatial distance traveled by the photon as measured in the instanta- 
neous local rest frame of the comoving observer (and shouldn't be 
confused with distance in any global frame). The "comoving" dis- 
tance parameter sm is the same for all photons with identical paths 
r(^) and depends on the fluid velocity u{x) they intercept, but not 
on their individual wavelengths. When Sm is used Eq. M2\ becomes 



dr dh dti diA ^ dh 

— I — + a(SM,A)—— 

dsM or dsM Ofi j oA 



where «(.?«, A) is defined as 
AdA 

a(sM,A:) = - —, 

and is related to Mihalas's a^isM, v] by 



X.-.^|)/. + ..(23) 



(24) 
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h 

By using Eqs. l |17| l, ( |20| (. and i2H . a(sM, /I) is easily evaluated, 



r ar 



as are the characteristic equations, 



dr 

-1— = J(M+Pl 
aSM 



dSM 



7(1 -y") 



(25) 

I, 

(26) 

(27) 
(28) 



These are to b e compared with Eqs. (3.4a), (3.4b) and (3.9) of 
iMihalaj (11980^. We have arrived at the original transport equation 
obtained by Mihalas ( 1980), characteristics and all, but do not pro- 
pose solving Eq. l l23l l However, if we were to now follow Miha- 
las, we would solve the single equation, ds^ = (hlX)d^, and ob- 
tain the integrated characteristics from l |16l > and (TtJi. By first find- 
ing the geodesies and secondly deciding on what variables to use 
on phase-space, we obtain the characteristics by substitution. For 
phase-space coordinates used in Eq. ( |23t , the characteristic curves 
had three non-vanishing tangent vectors Eqs. ( I26t . ( I27t . and i28t . 
The procedure we use is to choose as many coordinates on phase- 
space as possible that remain constant along geodesies (lifted to 
phase-space). For radial flows one constant coordinate is chosen 
from b and n e.g., 6, and the other two are <f itself and X. The 
transport equation is of the form l |18l l with characteristics having 
only two non-vanishing tangents, i.e., only two coordinates change 
along any characteristic, f(f) = ^ and A{^) [see Eq. J20H. Only one 
coordinate would change if we chose rather than A. However, 
the differencing procedure described in §6 requires the comoving 
energy (i.e., A for us) to be used as one of the phase-space coordi- 
nates. 

It is clear from the above discussion that if an afflne parameter 
such as ^ is used, the characteristic curves are completely deter- 
mined and don't have to be constructed as integral curves of their 
tangents. If a non-affine parameter is used (e.g., sm) all that must 
be done is to relate it to ^ (e.g., by inverting = - /(«(-«) • p)d^) 
alo ng the geodesic and subst itute. For a numerical example of 
see iBaron & Hauschildtl ( l20oi) . The affine parameter ^ is changed 
to optical depth by a similar substitution. A single reference 
wavelength such as /Isu = 5000 A is usually chosen, making 
rfTstd = ±;t-)sid'^'**'0 This choice greatly facilitates the generation 
of the spatial numerical grid. 



STATIC SPHERICALLY SYMMETRICAL 
SPACETIMES WITH STATIONARY FLOWS 



The relevant spacetime is 

ds^ = -e2*'''c2dr2 + e^'^^'^dr'- + r\de'- + sin^ ( 



(29) 



which we assume is asymptotically flat i.e., A(oo) = 0(00) = 0. We 
assume boundary conditions similar to the flat case, i.e., beyond 
r = R the fluid becomes transparent = x,i = 0) and that I(x, p) 
vanishes for incoming photons. With a stationary, but otherwise ar- 
bitrary fluid flow, the effective dimension of phase-space reduces 
to 6 (3 space and 3 momentum). For this case we also distort u{x) 



to coincide with static observers, i.e., (r, 0, 0)=constants, at some 
finite r value beyond r = R. 

The metric J29I ) has 4 Killing vectors, 1 time translation and 
3 rotations, which aid in finding photon orbits. In Fig. 2 we show 
the 1-parameter family of orbits confined to the 6 = n/2 plane and 
oriented symmetrically about ip = 0. They are labeled by the single 
impact parameter b and the wavelength /i„ seen by a rest observer 
at r = 00. The 3 non- vanishing momentum components are 



, dct 



n 



" g-2<t(r) 

h fee-*'"' 

dr_ 
d^ 



(30) 
(31) 
(32) 



y^2g-2[*(r)-<I>(i)] _ yi^ 



from which the spacetime orbits are 

^gA(r)+*(i)-2*(r) 



ct(r) 



6 y^2g-2[*(r)-*(i)] _ yi 

b e^<''> dr 



r be''" 



)-<bih)] _ 1,2 



(33) 
(34) 



All other photon orbits are obtained from these by rotations. Ac- 
cording to Euler, any active rotation may be described using three 
rotations R-((f)p + n/2)Ry(di,)R,(ifii,). For our purposes, (0p, 9p) spec- 
ifies the direction of the normal to the photon's orbital plane n,,ii„„, 
and ifib specifies the direction of the impact vector b within the orbit 
plane (see Fig. 2). These directions are space-like with no t com- 
ponents and are constrained by np/o,K • b = 0. Because the rotated 
orbit is orthogonal to n,,;„„p 



cos[0(r) - 0,,] = - cot dp cotS(r). 



The actual orbit in parametric form (r, 6(r), ifi(r)) is given by 



cos 9(r) 
tan 0(r) 



sin 9p s'm[(l>„,2(r) + 
cos Op tiin[tf>^,2(r) + 1//,,]- cot tf>p 
1 + cos 9p cot <pp tan[0;r/2(r) + if//,] ' 



(35) 

(36) 
(37) 



where (f>nii{r) is given by Eq. l l34t . Equation l l35t can be used in 
place of Eq. l l37t if desired. For an arbitrary stationary flow, the 
only net symmetry is a time translation. The specific intensity de- 
pends on six variables e.g., /( = /^(r, p). Instead of using the spher- 
ical polar angles and spherical polar momentum coordinates, we 
use the radius r, the comoving wavelength A, and four constants. It 
is convenient to choose the 4 constants from the set of 5 given by 
the impact vector b and the orbit plane normal e.g., the im- 

pact parameter b and the 3 angles ((f)p. Op, if/t,) described above. We 
now have = I^i(r,iipia„e,h), where the A dependence is implied. 
We could have eliminated A in favor of the constant Ac^, however, 
as with the flat-space case, resolution of the atomic lines and accu- 
rately calculating the angular integral [see Eq. l |47H dictates that we 
use the comoving wavelength. The transfer equation is still equa- 
tion l |18t , but now 



^ In classical plane-parallel and spherically symmetric radiative transfer r 
is measured from the outside inward and the minus sign is appropriate. 



dr_dh 

d^ dr 



(38) 
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and the single coordinate part of the characteristic to solve is 
Eq. ( I33t . The A{^) part is found by substituting into 



hydro calculation. For a numerical application of this section to the 
static Schwarzschild case see .Knop et aL (2007,) . 



y,.2g-2[*(r)-*(ft)J _ ^2 



If 



+ - sin Op sin[(^(r) - 0p]y3 



b cos 6,, 



r sin d{r) 



If 



(39) 



where the unit comoving 4-velocity u{x) has been written in terms 
of an otthonormal tetrad adapted to the static spherical polar coor- 
dinates of Eq. l |29l l. 



u(x) 



7 



d 



c dt 

1 d 
r sin 6 dd 



1 a 



(40) 



We assume that uix), and hence the /3'(r,d,(fi) are given by prior 
hydrodynamical calculations. Again as was illustrated Fig. [T] if /I 
is being used as a coordinate, it varies along a characteristic. 

Because of the paucity of hydro calculations in GR one can 
use yS's that come from Newtonian gravity calculations as a rea- 
sonable approximation as long as the gravity remains weak. The 
used in Eq. l l40b should be Newtonian components relative to a 
static spherical polar coordinate system. The GR metric associated 
with the Newtonian calculation is given by 



1/1- 



2Gm(r) A, 



(41) 



here Aq is the cosmological constant, M(r) is approximated as the 
Newtonian mass contained in radius r, and Pr(r) is the radial com- 
ponent of the pressure as seen by a static observer in the Newtonian 



6 LOGICAL STEPS FOR NUMERICAL INTEGRATION 

The solution of the spherically symmetric transfer equation for ra- 
dially moving flows has been discussed in detail by Hauschildt and 
Baron (Hauschildt 1992; Hauschildt &_Baron 1999, 2004b a). In 
particular iHauschildt & BaronI ( l2004al) showed how to stably dif- 
ference the dI,i/dA term. This method will also work in the more 
general case discussed here. We bri efly describe the app roach in 
simple terms, for more details see IHausch ildt & Baron (2004iJ). 
One first selects a fixed set of comoving wavelengths, A,„ at which 
to evaluate the specific intensity /„,, and treats dI,i/dA as a differ- 
ence, e.g.. 



dA 



(42) 



This turns the PDF, Eq. (TS}, into a discrete set of coupled 
DDEs with a single differential variable ^ for the set of /„,. The 
choice of the set [A,,,} is dictated by the variation of with A. Re- 
arranging, Eq. dlSI l becomes 



5A„, 



■ A„,-\ 



h A„ 
1'nj- + ] 



(43) 



(see'Mihalas"l980^ where A„, is determined by differentiating equa- 
tions l llSt or l |39l ) with respect to the affine parameter £. Even 
though phase-space is still 6-dimensional we are now attempting 
to find only on a discrete set of lines corresponding to constant 
comoving wavelengths A,,,. For both flat and static spherically sym- 
metric spacetimes, the remaining continuous variables are 3 posi- 
tion and 2 momentum coordinates. To simplify solving Eq. l |43l l, 
coordinates should be adapted to the particular spacetime symme- 
try and even then there is significant flexibility. For flat space we 
chose r and 2 of the n components making Eq. il6\ the characteris- 
tics. In Fig. 1 a spacetime geodesic of Eq. il6l is shown as it enters 



Steps for Solving the Radiative Transfer Equation for Arbitrary Flows in Stationary Spacetimes 7 



and exits the boundary r = R. The related characteristic curves in 
phase-space corresponding to a discrete set of wavelengths are also 
shown. If /i„ is used as the momentum coordinate, the character- 
istics of Eq. l llSt are all straight lines; however, if the comoving A is 
used the characteristics deviate from being straight but return some- 
where outside the boundaries where the comoving fluid coincides 
with the rest observers. The A part of the characteristics of the dif- 
ferenced transfer equation l l43t has also changed, i.e., comoving A„, 
now remains constant. The differencing term in Eq. I l42t is an ap- 
proximation which attempts to account for the change from /1„ = 
constant curves to /I = constant curves along which is propa- 
gated. In practice the differencin g procedure used is more compli- 
cated jHauschildt & Barorj|2004iJ) . When Eq. (TDl is solved for /,, 
at an exiting wavelength, e.g.,/l,„, its value depends on values 
along its prior path for a continuous spectrum of neighboring A val- 
ues. However, when Eq. J43l > is solved for /„, at an exiting point, 
its value depends on /,„/ values for only a discrete set of neighbor- 
ing wavelength A,„i. The discrete coupling appears in the RHS of 
Eq. l l43t in 77,,, and explicitly as 

The phase-space picture for the static spherically symmetrical 
spacetime is quite similar. We chose r and 4 constants from xipiane 
and b (see Fig. 2) and the only non-constant component of the char- 
acteristic curves of Eq. l l43t is given by Eq. l l33b . 

Often the source term is the more complicated part of the 
transfer equation. It contains various moments of the distribution 
function depending on the assumed physics of the photon-matter 
interactions. A detailed discussion of the moment f ormalism for 
steady state transfer can be found in iThornd jl98lh . To balance 
the interaction between the radiation field and the material, i.e., to 
obtain energy-momentum conservation of the radiation and mat- 
ter, one must adjust the material's parameters such as temperature. 
The opacity and Planck function depend on the temperature of the 
material and thus this temperature has to be consistent with the 
rate energy is being deposited by the photon gas (which is not in 
thermal equilibrium with the matter). One obtains this consistency 
through use of the energy-momentum tensor of the photon gas. Per 
unit comoving wavelength it is defined as the comoving solid angle 
integral 



rf(x) 



p) 



c/^.V^(., 

= ^ J"/,i(x,p)« + ;/")(;j(| + !/)dfl, (44) 

where we have decomposed the photon 4-momentum into two 
parts, one along the observer's 4-velocity a, and another perpen- 
dicular to it, u ■ 11,1 = 0, by defining 



p = -(p . u)(n„ + u). 



(45) 



Equation ( I44I I can be decomposed into energy, momentum, and 
pressure densities per unit wavelength as seen by u{x), 

Tf(x) = fEA(x)ii"(x)i/(x) + 



-f"^(x)i/(x) + -f,(x)u''(x) + pf(x), 
c c 



by defining 

r,(x) = j IAx,pKdn, 

p"f(x) ^ iA(x,pMdn, 



(46) 

(47) 
(48) 
(49) 



where = 0, pfu^ = and gapP^f = Q = '^^Ja/c, see Eq. l ll2b . 
The rate per unit comoving wavelength per unit volume that 4- 
momentum (energy/c and momentum) is being transferred to the 
photons is 



47T 

+ —B, 

c 



An 1 

— KABA-JAW--XAfA. 

c c 



(50) 



The familiar statement of radiative equilibrium (total absorptions 
equals total emissions in steady-state) for a comoving observer is 
the va nishing of the integral of u^T"^.^ over all wavelengths. Fol- 
lowing lLindauistl i 19661) it is convenient to define the particle (pho- 
ton) flux 4- vector 



r dP 



A_ 

ch 



eAx)u"{x) + -f^ix) 



(51) 



in terms of which the rate per unit comoving wavelength that the 
photon number density is changing due to absorption and emission 
by sources can be computed as 



An A 

— K,-[B,i-J,] 

c h 



(52) 



By inserting Eq. | |5U into the left hand side of Eq. J52b we obtain 
the identity 



(eAx)u'-{x)\, + -r,{x) 



An 



-K, [B, - 7,,] . 



(53) 



Integrating Eq. J53t over A, together with radiative equilibrium 
from Eq. l |50t , leads to the familiar result that the divergence of 
the flux is zero for a static fluid (i.e., for a fluid where u(x) oc the 
timelike Killing field). Equations l |50t and l |53l l are used to enforce 
the conditi on of radiative equilibrium (energy conservation, for ex- 
ample see iHauschildt & Barorj[l999t) . At depth, it is necessary to 
use the vanishing of the LHS of Eq. l |53t when integrated over A, 
rather than the RHS, because at high optical depth is very close 
to 6,1 and both are very large. Numerically, it is difficult to obtain 
an accurate result from the subtraction of two large numbers. 

The above expressions for the decomposition of T"^ and A'" 
are valid in any coordinate system and only require knowledge of 
Li(x) in that coordinate system. Choosing a comoving frame for an 
comoving observer in a non-symmetric spacetime essentially intro- 
duces an arbitrary rotation at every point in space. Fortunately it is 
not necessary to pick such a frame. When evaluating the comoving 
solid angle integrals in Equations l l47t -l|49b using stationary coor- 
dinates one eliminates the comoving element dCi. in favor of the 
stationary solid angle dClf, using 



dQ ■■ 



Up ■ p 
u ■ p 



dQn 



(54) 



Here mq is a unit vector pointing in the stationary frame's / direction. 
For the flat space case uq ■ p = -h/A„, u- pis given by Eq. l ll5t , and 



dQ 



-dfln 



For the spherically symmetric gravity field 



llQ 



(55) 



(56) 
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and M • p is given by Eq. l |39l l. To evaluate dO. we might have been 
tempted to introduce two comoving momentum variables yu, the 
cosine of the comoving polar angle, and oj, the comoving azimuthal 
angle, as customary. However, because of the arbitrariness of the 
flow evaluating these variables along a characteristic would have 
added two more complicated characteristics to solve and done little 
to aid integrating Eq. l l43t . 

Solving the transfer equation for an arbitrary stationary space- 
time is similar to the above flat space and spherically symmetric 
examples. A stationary spacetime is one that has a timelike Killing 
vector If the Killing vector is irrotational or equivalently hyper- 
surface orthogonal, the spacetime is called static as both examples 
were. Stationary coordinates (/, x') can always be adopted to the 
Killing vector, e.g., K = d/dt, which make the metric components 
independent of t. If the space is static, all dt®dx' cross terms of the 
metric can additionally be made to vanish by appropriately choos- 
ing a new t coordinate, t ^ t + fix'). When stationary coordinates 
are used, a fixed spatial boundary can be drawn around the source, 
beyond which photons no longer interact with the comoving fluid. 
The boundary is 2-dimensional and can be somewhat more compli- 
cated than the r = R of the two examples. Incoming photons can 
be labeled uniquely by 5 constants, 2 from the position they strike 
the boundary, 2 from the impact orientation with which they strike, 
and 1 from their wavelengths at oo (assuming the space is asymptot- 
ically flat). The exact position of any photon within the interaction 
region is then given by these 5 constants and the affine parameter ^ 
via the geodesic equations for the given geometry, and the comov- 
ing wavelength is given via Eq. {3}. Finding coordinates which are 
constants is also the essence of Hamilton- Jacobi theory of classical 
particle mechanics. A family of canonical transformation is sought 
on phase-space which takes the particles in phase-space from their 
initial positions to their positions at time t. The particles keep their 
initial values and only time changes. Solving Eq. {43} for the gen- 
eral case then proceeds exactly as in the examples. In the two ex- 
amples we used the existing symmetry of the spacetime to choose 
constants to label the photons rather than the boundary impact co- 
ordinates and angles. Impact constants could have been used; they 
can be computed from the symmetry constants we actually used, 
i.e., from b and n. However, it is easier to stick with symmetry con- 
stants when they exist. 



7 DISCUSSION 

Our formulation of the radiative transfer equation in terms of co- 
moving wavelengths and stationary coordinates, and the recogni- 
tion that the momentum directions can be pre-chosen by constants 
is the fundamental result of this paper. Schinder & Bludman (198?) 
recognized this for the case of purely static (no flow) transfer in 
spherical symmetry. Since the directions of the geodesies may be 
chosen for example at the boundary, the solution of the full 3-D ra- 
diative transfer problem in the presence of arbitrary hydrodynamic 
flows is very similar to the purely static case (no flow) described for 
example bv lHauschildt & BarorJ In that method space is di- 

vided up into a 3-D rectangular grid and long characteristics are fol- 
lowed from the outer boundary through the computational domain. 
The directions that the characteristics followed were simply chosen 
by dividing up the angular space at the boundary into equal parts. 
Since we have shown that the momentum directions may be chosen 
b y constants in both the flat-sp ace and curved-space, the procedure 
of lHauschildt & BaronI 12003) can be used with only the modifica- 
tion that Eqs. M5\ and J54l ( have to be evaluated once at each grid 



point. Naturally, detailed numerical tests have to be performed to 
ensure that the material properties such as density and composition 
are well resolved, and that there are enough "angles" to resolve the 
momentum-space variation of This differs from more classical 
methods in which the equations of the characteristics are solved in 
phase-space, with all momentum coordinates co-moving. 

If horizons are present in the spacetime, modification of the 
procedure outlined in this section may be necessary. However, for 
most astrophysical systems the emission of the radiation and the 
radiative transfer occur in the accretion disk (or in winds and jets) 
i.e., outside of any spacetime horizon. The obvious case where one 
would like to calculate both photon and neutrino transport in the 
presence of horizons, would be the format ion of a coUapsar, see 
for example iMacFadven & WooslevI h999t) : however, this would 
involve general relativistic MHD with radiative transfer, a feat that 
is far beyond current computational capabilities. 

Extending 1-D transfer calculations (e.g., spherical symmetry 
with radial flows) to 3-D applications with arbitrary flows is cur- 
rently of wide interest because of the desire to include the effects of 
rotation and Rayleigh-Taylor instabilities into stellar wind models, 
supernovae, and gamma-ray bursts; as well as the rapid improve- 
ment in computing capabilities which makes the extension possi- 
ble. Some confusion exists in which 1-D structures can be extended 
to 3-D and which cannot. We have pointed out that the wavelength 
A (or equivalently frequency) seen by a comoving observer u{x) is 
one that not only can be used, it is perhaps essential. Whereas a 
comoving frame is useful in the 1-D case, it should not be used 
in 3-D. Not that a comoving frame can't be defined for a comov- 
ing observer, it's just that there are too many possibilities and no 
natural way to make a unique choice between them. Consequently, 
as many as three additional functions of position (e.g., three angles 
at each point of spacetime) will be unnecessarily introduced into 
the transfer equation by introducing a comoving frame. The proce- 
dure we advocate for solving non-symmetric transfer problems is 
(1) start with the given spacetime geodesies, (2) change to appro- 
priate coordinates on phase-space, and (3) solve the transfer using 
Eq.(34). 

We have concentrated on stationary flows in stationary space- 
times (steady-state) applications, emphasizing the fact that, in 
spite of how coordinates might be chosen, characterist ics are re- 
ally geode sies extended to phase-space as described by iLindguislI 
fl966) and lEhlerj ( Il97ll) . We have also shown that by choosing an 
appropriate set of parameters (coordinates) to label the geodesies, 
characteristics can essentially become straight lines. When one of 
the coordinates is chosen as comoving wavelength A, only that par- 
ticular coordinate changes non-linearly. This slight complication is 
rewarded by a resulting simplification in the absorption and emis- 
sion terms on the "collision" side of the transfer equation. Addi- 
tionally the use of A allows one to convert the transfer equation 
from an integro-PDE to a system of integro-ODE's by a differenc- 
ing procedure which restricts comoving wavelengths to a discrete 
set {/!„,). These wavelengths remain constant along the characteris- 
tics of the now differenced transfer equation and greatly simplifies 
constructing the formal solution using long or short characteristic 
methods. 

Other authors have recognized the current need for adapting 
1-D m ethods of solution to 3-D problems in radiative transport. 
ICarda ll et al. (2005) derive the relativistic equation of transfer in 
flat spacetime, with similar goals to this paper; however, their ap- 
proach is somewhat geared to the discrete ordinates matrix expo- 
nential method of numerical solution, or to generalized variable 
Eddington factor (moment) methods. 
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iBrodericS j2005h also recognized the advantage of using con- 
stant momentum variables, which he attempts to define through 
the use of Fermi-Walker transported tetrads. Even though his use 
of these tetrads isn't clear, parallel transport itself results in con- 
stant momentum components. If the spacetime is curved, parallel 
transporting along every geodesic results at an infinite number of 
tetrads at every point. And, if scattering exists, the relation between 
these tetrads must be ascertained before the emissivity integral [see 
Eqs. jl2b and J47H can be evaluated. 

In conclusion, we have presented a workable outline of solv- 
ing the radiative transport equation for many 3-D steady-state prob- 
lems. 
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